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Abstract: The purpose of this paper is to present a simple micromechanics-based model to 
estimate the effective thermal conductivity of macroscopically isotropic materials of matrix- 
inclusion type. The methodology is based on the well-established Mori-Tanaka method for 
composite media reinforced with ellipsoidal inclusions, extended to account for imperfect 
thermal contact at the matrix-inclusion interface, random orientation of particles and parti- 
cle size distribution. Using simple ensemble averaging arguments, we show that the Mori- 
Tanaka relations are still applicable for these complex systems, provided that the inclusion 
conductivity is appropriately modified. Such conclusion is supported by the verification of 
the model against a detailed finite-element study as well as its validation against experimen- 
tal data for a wide range of engineering material systems. 

Keywords: Engineering materials; homogenization; Mori-Tanaka method; thermal conduc- 
tivity; imperfect interface; finite element simulations 



1. Introduction 

There has been a clear trend over the last decade to exploit ever greater detail of the material structure 
towards better predictions of its response from simulations. Hierarchical modeling strategy, regardless 
whether coupled or uncoupled but mostly of the bottom-up type, has served to provide estimates of the 
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macroscopic response. In this process, geometric details decisive for a given scale are first quantified 
employing various statistical descriptors [1], but eventually smeared via homogenization to render larger 
scale property. Greater precision is expected when introducing the results of microstructure evaluation 
into the homogenization step. However, the actual gain when compared to the cost of this analysis is still 
in question. Obviously, description of evolving microstructures or rigorous representation of deformation 
mechanisms would require to account for almost every detail of the microstructure on a given scale. But 
how deep do we have to go if only the effective macroscopic response (i.e. linear macroscopic properties) 
is of the primary interest? Such a goal is addressed in this contribution. 

Here, the modeling effort concentrates on the evaluation of effective thermal conductivities of var- 
ious engineering materials with a significant degree of heterogeneity whereas focusing on imperfect 
thermal contact along constituents interfaces. We shall argue, shielded by available experimental data, 
that reasonably accurate predictions of macroscopic response can be obtained with very limited infor- 
mation about actual microstructure such as volume fractions and local properties of material phases. 
Consequently, we lump the entire analysis on the assumption of representing true material structures by 
statistically isotropic distribution of spheres. Figure 1 shows micro-images of selected material repre- 
sentatives which seem to admit this classification. Note that whatever material phase embedded into the 
matrix (the basic material) is henceforth termed the heterogeneity in real material systems while it is 
termed inclusion in approximations adopted for calculations. 

Figure 1. Examples of micro-graphs of real engineering materials taken in back scattered 
electrons: a) Alkali- activated fly ash, b) Alumino- silicate ceramics with Fe and silicium par- 
ticles (dark phase), c) Superspeed - alloying ingredient into crude iron for cast iron working 
with silicon particles (dark phase). Reproduced with permission of L. Kopecky (CTU in 
Prague). 




Strong motivation for this seemingly swingeing simplification is supported by experimental measure- 
ments presented in [2] for cement matrix based mixture of rubber particles and air voids. Comparison 
between experimental data and predictions provided by the Mori-Tanaka averaging scheme under the 
premise of random distribution of spherical inclusions, the method in this particular format promoted 
herein, appears in Figure 2. The match is almost remarkable. 

Going back to Figure 1 one may object that while admissible for Figure 1(a) the spherical approxi- 
mation of particle phase in Figure 1(c) will yield rather erroneous predictions. Note, however, that this 
attempt is not hopeless providing the microstructure can still be considered as macroscopically isotropic, 
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Figure 2. a) Evolution of effective thermal conductivity x as a function of volume fraction 
of rubber in solid phase, b) correlation of measured and calculated values; p is the correlation 
coefficient. 
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ensured for statistically isotropic distribution of heterogeneities having isotropic material symmetry. In 
that case, it can be shown that the previously mentioned Mori-Tanaka method written out for spherical 
inclusions is adequate providing the material properties of the inclusions are suitably modified. Al- 
though this step requires information beyond that of volume fractions of phases, the benefit of gathering 
additional data will become particularly appreciable once turning our attention to material systems with 
imperfect interfaces, the principal objective of this study. 

The problem of quantifying the influence of imperfect thermal contact on the overall thermal conduc- 
tivity has been under intense study in the past. Hasselman and Johnson [3] provided estimates for dilute 
concentration of mono-disperse spherical and cylindrical heterogeneities. Successful application of this 
simple model to Al/SiC porous composites is presented in [4]. The Hasselman- Johnson results were 
then extended by Benveniste and Miloh [5] to spheroidal particle shapes with imperfect interfaces and 
subsequently applied in the framework of the Mori-Tanaka method [6]. These early developments were 
later generalized by Nogales and Bohm, who proposed in [7] a simple method for dealing with polydis- 
perse systems of spherical particles. In addition, rigorous third-order bounds for effective conductivity 
of macroscopically isotropic distribution of particles with imperfect interfaces were derived by Torquato 
and Rintoul [8]. Alternatively, as demonstrated by Hashin [9], the material systems with imperfect inter- 
faces can be accurately approximated by the coated inclusion model due to Dunn and Taya [10], which 
also accounts for different orientation of the inclusions. If limiting attention to spherical inclusions the 
results presented in [7] can be obtained in a very elegant way by simple extension of one-dimensional 
analysis. This is demonstrated in Appendix B. 

To exploit this result in practical applications of the Mori-Tanaka method to a heat conduction prob- 
lem, prediction of effective thermal conductivity in particular, we adopt the analysis scheme graphically 
presented in Figure 3. We start from the assumption of multidisperse system of randomly oriented 
spheroidal inclusions with possibly imperfect thermal contact (non-zero temperature jump along the in- 
terface). To arrive at the desired approximation of multidisperse system of spherical inclusions with 
perfect interfaces (temperature continuity along the interface) we proceed in five consecutive steps. 
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Figure 3. Mori-Tanaka based scheme: Strategy of derivation 




(iii) Orientation averaging (iv) Imperfect interface (v) Polydispersivity 




These steps are mathematically described in Section 2. Section 3 is devoted to both validation and 
verification of the proposed scheme against available experimental data and finite element simulations 
performed for several representatives of statistically isotropic random microstructures. The crucial re- 
sults and principal recommendations are finally summarized in Section 4. 

2. Theoretical background of the Mori-Tanaka method 

In this section attention is accorded to essential theoretical details of the Mori-Tanaka method in view 
of the five steps in Figure 3. In the first step, we consider a single inclusion with perfect interface sub- 
ject to far- field loading. This step is theoretically elaborated in Section 2.1. Solution of this problem is 
then employed in Section 2.2 to estimate the overall conductivity of a composite consisting of multiple 
ellipsoidal inclusions bonded to a matrix phase. The third step addressed in Section 2.3 is reserved for 
systems with randomly oriented inclusions with a uniform distribution over the hemisphere. Here, a sim- 
ple orientation averaging argument is shown to demonstrate that the effective conductivity of the system 
coincides with the conductivity of a system reinforced by spherical inclusions, thermal conductivity of 
which is appropriately modified. Following [7], an analogous argument is employed in the next step 
outlined in Section 2.4 to account for imperfect thermal contact along the matrix-inclusion interface. It 
is shown that in this case the modified conductivity becomes size-dependent. This eventually allows us 
to extend the scheme to polydisperse systems in Section 2.5. 

2.1. Single inclusion with perfect interface 

Let us consider an ellipsoidal inclusion Vt l , with semi-axes a\ < a 2 < a 3 embedded in the matrix 
domain f2 m . We attach to the inclusion a Cartesian coordinate system with the origin at the inclusion 
center and axes aligned with the semi-axes. The distribution of local fields follows from the problem 



q(x) = -x(x)h(x), V T q(x) = fora?eIR 3 , 



(1) 
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where qeM 3 denotes the heat flux, h G M 3 denotes gradient of the temperature 6 (i.e. h(x) = *V9(x)) 
and x designates the 3x3 symmetric positive-definite matrix of thermal conductivity given by 

X\ x ) = \ (2) 
I x m otherwise. 

Eq. (1) is completed by the far- field boundary conditions, cf. [5, Eq. (14)], 

9(x)=H T x for ||a;|| -> oo, (3) 

with H G IR 3 denoting the overall (macroscopic) temperature gradient. Due to linearity of the problem, 
we can introduce the temperature gradient concentration factor A G IR 3x3 in the form: 

h(x) = A(x)H for x G K 3 . (4) 

As shown first by Hatta and Taya [11], the concentration factor is constant inside the inclusion and admits 
the expression 

A-\x) = {A')' 1 = I-S (x™)- 1 (X m - X 1 ) for x G Q\ (5) 

where I denotes the unit matrix and S G IR 3x3 is the Eshelby-like tensor which depends only on the 
matrix conductivity x m an d the ratios of semi-axes lengths fa = a<i : a\ and fa = : a%, see also 
Appendix A for additional details. For the spherical inclusion with isotropic conductivity x 1 — X 1 ! 
embedded in the isotropic matrix with x m — X m Ii Eq. (5) simplifies as 



2x m + x i- 



A^A^I with Al ph = - S. r - (6) 



2.2. Multiple inclusions with perfect interface 



In the next step, we adopt the results of the previous section to estimate the overall behavior of a 
composite material consisting of distinct phases r = 0, 1, . . . , N. The value r = is reserved for the 
matrix phase f2 m and the r-th phase (r > 0) corresponds to the ellipsoidal inclusion flM, characterized 
by its semi-axes o!(' ,c4 and Og , volume fraction and conductivity x^ ■ Following Benveniste's 
reformulation [12] of the original Mori-Tanaka scheme [13], the interaction among phases is approxi- 
mated by subjecting each inclusion separately to the mean temperature gradient in the matrix phase H m 
in Eq. (3). As a result, the temperature gradient inside the r-th phase remains constant and reads 

H (r) = T (r) H m fQJ . r = 0, 1, . . . , jV. (7) 

Here, 

and is the 3x3 partial temperature gradient concentration factor of the r-th phase, given by 

/• J * forr = 0, 

" R^T { [\R^) T forr = l,2,...,iV, { ) 
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where the 3x3 rotation matrix accounts for the difference in the global and local coordinate 
systems, see Section 2.3 for additional details, and Tp equals to A 1 in Eq. (5) with x 1 — an d S 
determined from values a(' , c4 and a^. 

The volume consistency of the overall temperature gradient H and local averages requires 

N / N \ 

H = J2 c (r)H(r) = c (r) T (r) I H m . (10) 

r=0 \ r=0 / 

Inverting the above equation gives the average temperature gradient in the matrix as 

H m =[^c (r ¥ r) | H = A m H, (11) 



r=0 



which, when substituted into Eq. (7), yields the explicit expression for the phase temperature gradients 
in the form 



A' 



H (r) = T (r) ^ c (r) T (r) jj = A (r) H ^ (12 ) 



v r=0 



where A m , A^ r ' are the matrix and inclusion temperature gradient concentration factors, respectively. 
As each phase is assumed to be homogeneous, the average heat flux in the r-th phase 

Q (r) = ToTtY! I dx for r = 0, 1, . . . , N, (13) 

equals to 

q(0 = _ x (r) H (r) for r = 0, 1, . . . , iV. (14) 
This allows us to express the macroscopic heat flux in the form 

N / N \ / N \ _1 

q = j2 c (r) Q (r) = - E c (r) x (r) T« Yl c(r)T(r) H w 

r=0 \r=0 / \r=0 / 

from which we obtain the effective conductivity Q = —x H H in its final form 

X R = (^c m x m + c (r) x (r) T (r) j |Vl + c(r)TW j • ( 16 ) 

Finally, assuming that the composite consists of isotropic matrix with conductivity x m an d spherical 
isotropic inclusions with conductivities x^ r \ Eq. (16) becomes x R — X U I> with 

c m X m + f> (r) X (r) T£ 

< * r Q III 

* H = '-T - d ^"amj- < 17 > 

r=l 
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2.3. Orientation averaging 



We are now in a position to provide estimates of the effective thermal conductivity for composites 
with M (with M <C N) inclusion classes indexed by s = 1,2,..., M. Each class is characterized by 
a single Eshelby-like matrix S in Eq. (5) and represents the reference ellipsoidal inclusion randomly 
oriented over the unit hemisphere with an independent uniform distribution of orientation angles. 

To this goal, consider a quantity X e £ IR 3x3 , expressed in a local coordinate system aligned with a 
certain reference inclusion. Its form in the global coordinate system follows from 



where cp, <p and ip denote the Euler angles 1 and the transformation matrix R is provided by 



(18) 



R{y,<p,ip) 



cos ip — sin ip 
sin ip cos ip 
1 



cos 4> — sin <p 
1 
sin <b cos <t> 



cos ip — sin (p 
sin ip cos 
1 



(19) 



The orientation average of X is denoted by double angular brackets: 

"27T /*7T /*27T 







^0 



X(ip,(j), ip) sin d0 dip. 



(20) 



Straightforward calculation, presented e.g. in [14, Appendix A. 2. 3], reveals that the orientation averag- 
ing of an arbitrary X? £ IR 3x3 yields 



((X)) = ((X))I with «X)) = ^PQ 



(21) 



8=1 



Repeating the steps of the previous section with partial temperature gradient concentration factors 
replaced with their orientation averages, we obtain, after some manipulations presented e.g. in [15, 
Section B], the scalar homogenized conductivity in the form 



M 



c m X m + $>M(( X ( s )T 



(*)\ 



X 



s=l 



M 



(22) 



c m + ^c (s) ((T (s) » 



s=l 



Therefore, it follows from the comparison of Eq. (22) with Eq. (17) that the system of randomly oriented 
inclusions embedded in an isotropic matrix is indistinguishable, from the point of view of homogenized 
conductivity, from the system of spherical inclusions with an apparent conductivity 



X 



Q) _ ix {s) T^)) 
«7»}} ' 



(23) 



'Note that so-called "x2 convention" is used, in which a conversion into a new coordinates system follows three consecu- 
tive steps. First, the rotation of angle ip around the original X 3 axis is done. Then, the rotation of angle <f> around the new x 2 
axis is followed by the rotation of angle tp around the new X3 axis to finish the conversion. 
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which yields 

M 

with &l = 

2.4. Imperfect interface 



y H = !=* with T (s > = — (24) 



sph 

s=l 



The presence of imperfect thermal contact at the matrix-inclusion interface dVl 1 results in temperature 
jump [0], the magnitude of which is provided by Newton's law, e.g. [16, Section 1.3]: 

n T (x)q(x) = k(x)lO(x)j forxedW, (25) 

where k denotes the interfacial conductance (with k — > oo corresponding to perfect interface and k = 
to ideal insulation) and n denotes the normal vector at the interface oriented outside the inclusion. This 
relation, together with Eqs. (l)-(3), defines the single inclusion problem accounting for the presence of 
imperfect interface. Its solution is, however, substantially more involved as the temperature gradient in- 
side an ellipsoidal inclusion becomes position-dependent; the concentration factor is then available only 
in the form of complicated infinite series expansion for spheroidal inclusions [5] or ellipsoidal coated 
inclusions [10]. Nevertheless, when restricting the attention to spherical inclusion of radius a, it can be 
shown that the temperature gradient within inclusion recovers the constant value and the concentration 
factor becomes [5] 

-~ 3v m ■ ■ ak 

A s P h(x\k,a) = ————— — r where xKx\ k, a) = x'-j—, — r, (26) 

2x m + X 1 (X 1 ,M) ak + x 1 

see also Appendix B for simple derivation of this result. Hence, analogously to the previous section, the 

interfacial effect can be modeled by replacing the "true" conductivity \ l by the size-dependent apparent 

value x 1 provided by Eq. (26) 2 . Assuming in addition that each class of inclusions is characterized by 

identical semi-axes lengths and a!j and interfacial conductance k^ s \ we propose to extend the 

relation (24) into the form, cf. [7, Section 2] 

M 

c m X m + y c (s) X {s) T^l 

X M ' sph 9 v m i Cfs)' X X n(s)U(s) i yAd '> 



s=l 



with and the apparent conductivity x^ given by Eq. (23). 

2.5. Polydisperse systems 

Even though Eq. (27) is applicable to very general material systems, in practice we typically assume 
single inclusion family, with the particle size distribution characterized by a probability density function 
p(a) satisfying 

/oo 
p(a) da = 1. (28) 
-oo 
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In this context, the effective conductivity finally becomes 



c m x m + c (l) I £ (1) T (1) 

x = f^n ' (29) 

where, for an arbitrary function g(a), {g} denotes its expected value given by 

POO 

{g} = / g(a)p{a) da. (30) 



Following e.g. [7,17], the log-normal distribution with the probability density function 



p(a) = — exp 

V27racr 



ln(a) — \i 
V2a 



2 N 



a>0, (31) 



will be employed to characterize materials' polydispersity. The parameters ji and a are provided by [7, 
Eq. (17)] 



i / \ 1 i f S + VS 2 + A\ a 90 -a 10 
M = ln(a 50 ), - = ^ln i , S = , (32) 



where a x denotes the x-th percentile of the particle radii and S is the span of the size distribution. 
3. Mori-Tanaka estimates - example results 

3.1. Validation against available experimental data 

Two particular examples of real engineering materials are examined in this section to show applica- 
bility of Eq. (27) and its extension for polydisperse distribution of heterogeneities, Eq. (29), even when 
disregarding their actual shape and simply accepting a spherical representation of the inclusions in the 
Mori-Tanaka predictions. The results provided by these two equations are corroborated by available 
experimental data. 

Random dispersion of copper particles in the Epoxy matrix 

In this first example we compare the single-phase Mori-Tanaka predictions with the experimental 
results of de Araujo and Rosenberg [18] who measured the effective thermal conductivities of systems 
consisting of random dispersion of metal particles in the epoxy matrix for several values of the interfacial 
resistance arising due to acoustic mismatch at the particle-matrix interface particularly for low temper- 
atures below 20K. To enhance credibility of the Mori-Tanaka predictions we focus on one particular 
system made from epoxy resin filled with copper particles also examined in [8] in view of the three- 
point lower bound assuming a random array of superconducting hard spheres (i.e. x ~~ >* °°)- It has 
been shown, see [18,19], that for metal-filled composites with ratio a = x /x m > 10 2 the macroscopic 
conductivity does not depend on the thermal conductivity of the metallic filler, but only on its volume 
fraction together with the properties of matrix and matrix-particle interface. This becomes evident when 
rewriting Eq. (27)i in terms of a 

_ i _,_ 3cgj [a-(l + iZ)] 

X m c™[a + 2(l + R)}+3cW(l + RY K ) 
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where we introduced the dimensionless quantity R adopted in [8] 

R=\r = ^, (34) 

ka ka 

where a is the sphere radius. It follows from Eq. (33) that for the inclusion size equal to the Kapitza 
radius [20] 

m m 

« X _ X ~« 

ok = --r- « — (35) 

a — 1 k k 

the effective conductivity equals to that of the matrix, thus the effect of inclusions becomes completely 
shielded by the interface. For a < ax, the overall properties are dominated by interfaces and the effective 
conductivity decreases with increasing even when a > 1. For a > a K , the bulk properties of phases 
become dominant; see also [20] for further discussion. 

In addition to experimental measurements, we also present a comparison with the Torquato-Rintoul 
three-point lower bound in resistance case derived in [8, Eqs. (8)], evaluated for the statistical parameter 
C(c^) given for the model of impenetrable spheres in [21, Table II (simulations)]. Note that the inter- 
facial properties can be estimated by measuring the ratio of temperature jump to the applied heat flux 
across a thin bi-material layer, by acoustic mismatch model [4] or, indirectly, from an inverse approach 
as discussed below. The results are presented for two different temperatures 9 = 4K (a = 14.8ok) and 
= 3K (a = 4.93a K ). 

Figure 4. Evolution of the normalized effective thermal conductivity x R /x m as a function 
of a) phase contrast a and relative particle radius a/a K (c^ = 30%) and b) volume fraction 
of copper particles, c) correlation of measured and calculated values; p is the correlation 
coefficient. 




The influence of parameter a together with expected particle size dependence, now hidden in param- 
eter R through Eq. (34)i, is evident from Figure 4(a) plotted for 9 = 4K (a K = 3.38 /mi). Clearly, the 
Mori-Tanaka predictions confirm negligible influence of observed for particulate composites with 
X^ 3> X m ( a > 10 2 ) as well as decreasing trend for x H with decreasing particle size caused by imper- 
fect thermal contact. Figure 4(b) then displays evolution of normalized effective thermal conductivity 
as a function of volume fraction of copper particles. Predictive capability of the Mori-Tanaka method 
is supported here by a very good agreement with both experimental measurements and the Torquato- 
Rintoul bounds [8]. Finally, Figure 4(c) shows correlation between theoretical (MT predictions) and 
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experimental results. The solid line was derived by linear regression of measured and calculated ef- 
fective conductivities using the least square method. Another indication of the quality of numerical 
predictions can be presented in the form of Pearson's correlation coefficient written as 

n(xy) - (x)(y) 

P — I I ■> (3o) 

A/n(x 2 ) - (x) 2 y/n(y 2 ) - (y) 2 

where (a) = Ym=i a " n * s me num ber of measurements, x stands for the experimental and y for the 
corresponding theoretical values. Note that for p = 1 the correspondence is exact. In this case the 
correlation coefficient equals to 0.999 suggesting almost perfect match between measured and predicted 
values, also evident from graphical presentation. 

Al/SiC composite 

In [4] the authors studied the effect of imperfect thermal contact on the macroscopic response of 
Al/SiC porous composites. The paper presents the results of a thorough experimental investigation and 
traces of an inverse approach in material mechanics for inferring material properties of unknown com- 
ponents of the composite by matching numerical and experimental results. This approach was first 
exploited to derive the matrix thermal conductivity from known electrical conductivity of the composite. 
Next, the Hasselman and Johnson model [3] was employed under the premise of random distribution 
of spherical particles of identical size to estimate the particle thermal conductivity and interfacial ther- 
mal conductance for pore-free specimens and subsequently utilized in the two-step application of the 
Hasselman- Johnson model to address the influence of pores. Note that the material data used in these 
predictions can be thought as optimal with respect to the adopted Hasselman- Johnson model. 

Hereinafter, we compare the results presented in [4] for seven specimens with pore-free matrix. In 
addition, we take advantage of available characteristics of the SiC particles, the span S and the 50-th 
percentile a 50 , to extend the analysis to polydisperse systems as presented in Section 2.5. The input 
material data are listed in Table 1 . 

Table 1. Material properties [4]. 



Al matrix 


SiC particles 


Interface 


[Wm-'K" 1 ] 


[Wm^K" 1 ] 


187 


252.5 


72.5xl0 6 



The Mori-Tanaka predictions stored in Table 2 were provided by Eq. (29). The integral (30) was 
evaluated such that the entire interval was split into 1,000 segments, thus considering 1,000 different 
particle sizes of spherical shape. Within each segment s the probability function p(a) was approximated 
by a straight line and the volume fraction cW of a given set of particles, given by the segment area, 
was then associated in a logarithmic way with the mean radius of this set of particles. Standard 
trapezoidal integration rule was then used to sum over all 1,000 segments. Examples of probability 
distribution functions for specimens No. 1 and 7 are plotted in Figure 5. 

Graphical representation of the results is further seen in Figure 5(b),(c) confirming the sensitivity 
of the effective thermal conductivity to the mean particle size distribution. Note that individual points 
in Figure 5(b) correspond to slightly different volume fraction, see Table 2. 
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Table 2. Characteristics of SiC particles (cf. [4, Table 1]), and comparison of effective 
thermal conductivity x H [Win -1 K -1 ] between experimental measurements [4] and MT pre- 
dictions Eq. (29). 



Sample Radius \pm] SiC Results 



No. 


aio 


a 90 


S 


vol. 


Exp. 


MT 


1 


55 


114.5 


0.71 


0.58 


219 


217.8 


2 


23 


65.5 


1.02 


0.58 


210 


212.3 


3 


19.5 


37.5 


0.66 


0.60 


208 


208.5 


4 


11.5 


25 


0.79 


0.59 


198 


199.9 


5 


7 


17 


0.86 


0.58 


195 


190.8 


6 


5 


12 


0.82 


0.55 


184 


182.5 


7 


2.4 


7 


1.05 


0.53 


160 


161.3 



Figure 5. a) Examples of probability distribution functions of particle radii a for specimens 
No. 1 and 7, b) evolution of effective thermal conductivity x H as a function of particle radius 
a, c) correlation of measured and calculated values; p denotes the correlation coefficient. 




Almost negligible deviation from experimental results measured as 

£(x exp - x MT ) 2 

E = = 1.1% 

* £(x MT ) 2 

is, however, not surprising owing to the used material parameters, which were not measured but rather 
fitted to a micromechanical model. Even when comparing the Pearson correlation coefficients, 0.98 for 
the Hasselman- Johnson model and 0.993 for polydisperse MT model, the improvement when accounting 
for more accurate representation of particle size distribution is marginal. This can be explained by a 
very small variance associated with adopted distributions, recall Figure 5(a). Nevertheless, it is fair to 
point out that unlike the Hasselman- Johnson model the Mori-Tanaka approach is not limited to spherical 
particles providing the transformations given by Eqs. (24) and (27) are admissible. The influence of 
shape of particles on the macroscopic response has been put forward, e.g. by Jackel in [19], and is 
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numerically investigated in the next section suggesting increasing thermal conductivity of the composite 
with transition from spherical to needle-like particles, the trend also observed experimentally [18,19]. 

3.2. Verification against finite element simulations 

It has been argued in the previous sections that even very limited information about micro structure 
amounted to phase properties and corresponding volume fractions might be sufficient to provide a rea- 
sonable estimate of macroscopic response of various engineering material systems generally classified 
as being macroscopically isotropic. This naturally invites the assumption of spherical representation of 
otherwise irregular heterogeneities. Although supported by several practical examples discussed in the 
previous section, we should expect and even identify, at least qualitatively, limitations to such perception. 
In doing so, this section presents numerical investigation of some specific issues such as the influence 
of shape and size of inclusions or mismatch of phase material properties on the predicted macroscopic 
response. 

All numerical results reported bellow are obtained in the two-dimensional setting, hence the argu- 
ments presented in Section 2 need to be translated to the planar case. In particular, inclusions are mod- 
eled as elliptic cylinders of aspect ratio fi% = a 2 /ai with a 3 — > oo. The corresponding Eshelby-like 
tensor is given by Eq. (42), which leads to the two-dimensional thermal concentration factor for circular 
inclusion in the form 



A 1 



^cir / (2x2) With A' ch 



2x r 



The apparent conductivity due to random orientation is provided by 



(37) 



X 



M _ <(x w r (,) »gp 
«rw»2D ' 



(38) 



where ((•^d stands for the orientation averaging for the uniform distribution of the Euler angle (j). 
Finally note that the apparent conductivity due to imperfect interface is provided by the same relation as 
in the three-dimensional setting, compare Eq. (26) with Eq. (50). 

Figure 6. Examples of random macroscopically isotropic microstructures: a) circular cylin- 
ders (/3 2 = 1), b) elliptical cylinders with aspect ratio (3 2 = 3, c) elliptical cylinders with 
aspect ratio /3 2 = 9. 
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Three particular representatives, generated such as to approximately resemble the real microstruc- 
tures in Figure 1, appear in Figure 6. To comply with general assumptions put forward in the previous 
sections, we consider locally isotropic phases with variable contrast of material properties. Additionally, 
we assume the above microstructures being periodic and adopt the classical first-order homogenization 
strategy, see e.g. [22,23], to provide estimates of the macroscopic response. The results are plotted 
in Figure 7. 

Figure 7. Variation of the normalized effective conductivity x U /x m f° r three microstructures 
in Figure 6 with perfect interfaces, determined by periodic FEM homogenization for phase 
contrast a) a = 3, b) a = 10 and c) a = 20; /?2 denotes the inclusion aspect ratio and 
the inclusion volume fraction. 




These results clearly indicate not only the influence of the shape of inclusions on the macroscopic 
response but also a strong dependence of these predictions on the contrast of material properties of indi- 
vidual phases. Thus drawing from the plots presented in Figure 7(a) one may suggest that the proposed 
circular representation of generally non-circular heterogeneities is still acceptable when their shapes only 
moderately deviate from a circle and when the mismatch of phase properties is not too severe, which cer- 
tainly is the case of a number of real materials as demonstrated in the previous section. This expectation 
is quite important particularly when dealing with imperfect thermal contact in which case only spherical 
and circular inclusions can be easily handled analytically. 

If the circular approximation of heterogeneities is no longer acceptable or the contrast of phase prop- 
erties is excessive, one needs to look for more details about microstructure. In such a case, even two- 
dimensional images of real systems, at present almost standard input for any material based analysis, 
may play an important role in assessing better approximations of shapes, say elliptical, of these hetero- 
geneities, see e.g. [24]. Then, being given the elliptical shape of the inclusion allows us to appropriately 
modify its material data, recall Eq. (38), and define a certain indicator of the real microstructure k con , 
e.g. as a ratio of the modified and original partial temperature gradient concentration factors 

f (1) 

^'corr (39) 
"^cir 

now determined for circular inclusions. 
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Figure 8. a) Evolution of correction factor k corr for systems with perfect interfaces as a func- 
tion of the ratio of semi-axes (3 2 and variation of normalized effective conductivity x H /x m 
obtained by FEM and MT with modified conductivity x> l > for phase contrast b) a = 3 and 
c) a = 20. 
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Variation of this parameter as a function of the shape of inclusion is seen in Figure 8(a) further con- 
firming quite strong influence of the phase properties mismatch. 2 The modified conductivity when intro- 
duced successively into Eq. (37) 2 and Eq. (24)i then renders the estimate of effective conductivity almost 
identical to actual micro structure with non-circular inclusions as evident from plots in Figure 8(b),(c). 
Note that only the first and the third microstructure in Figure 6 were examined to first confirm that the 
Mori-Tanaka method is indeed well suited for statistically isotropic random microstructures and second 
to promote applicability of this simple transformation from elliptical to circular representations even for 
shapes markedly distinct from circles. Small but evident deviation of the results observed in Figure 6 
for the third needle-like microstructure and large mismatch of conductivities of the inclusion and ma- 
trix equal to 20 can be attributed to the finite size, although infinite in the sense of periodicity, of the 
representative model not large enough to yield statistically isotropic microstructure. 

Figure 9. Statistically isotropic distribution of circular cylinders with variable radius of their 
cross-section. 

(a) (b) (c) 



2 Note that the parameter fc corr is determined for two dimensional systems and thus is not applicable to ellipsoidal inclusions 
of the same aspect ratio. 
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The second set of numerical simulations addresses the theoretically derived dependence of macro- 
scopic predictions on the size of inclusions in cases with imperfect interfaces generating jumps in the 
local temperature field. For simplicity, we limit our attention to statistically isotropic distributions of 
circular cylinders with a radius varying from sample to sample. Three such microstructures are shown 
in Figure 9. Note that the same volume fraction and the same number of inclusions was maintained in all 
simulations. Zero thickness interface elements were introduced to account for imperfect thermal contact. 

Figure 10. a) Variation of normalized effective conductivity x H /x m f° r a system with im- 
perfect interface (k = 60 x 10 6 Wm^K- 1 and = 100 Wm^K -1 , x m = 10 Wm^K" 1 
and c- 1 ' = 40%) as a function of inclusion radius a and b) variation of effective conductivity 
for systems weakened by cylindrical voids as a function of inclusion volume fraction c^; (3 2 
denotes inclusion aspect ratio. 




The relevant results appear in Figure 10(a). Both the results found from finite element simulations and 
corresponding Mori-Tanaka predictions are displayed to clearly identify the mentioned size dependence. 
Proper modifications in the sense of Eq. (38) now become even more important as indicated by the 
results generated for elliptical microstructures from Figure 6 with cross-sectional area equal to the area 
of the circle. These are indicated by x-symbol and the ratio of semi-axes of the corresponding elliptical 
cylinder. The Mori-Tanaka estimates found from the application of Eqs. (38) and (26) are reasonably 
close further supporting the proposed approach for the modeling of real materials with an imperfect 
thermal contact. Intuitively, it can be expected that the value of interfacial conductance k may also 
show some effect as to the estimates of effective conductivities for non-circular inclusions. This notion 
is supported by the results presented in Figure 10(b) showing variation of effective conductivity of an 
isotropic matrix weakened by voids with a very low conductivity. Clearly, the influence of shape of the 
inclusions is quite pronounced. 

4. Conclusions 

The Mori-Tanaka micromechanical model has been often the primary choice among engineers to 
provide quick estimates of the macroscopic response of generally random composites. Motivated by an 
early theoretical as well as experimental works on this subject, the Mori-Tanaka method was examined 
here in the light of the solution of a linear steady state heat conduction problem allowing us to estimate 
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the effective thermal conductivity of variety of real engineering materials experiencing an imperfect 
thermal contact along the constituents interfaces. 

Adhering to the only limitation, an assumption of macroscopically isotropic composite, it was shown 
that the method originally proposed by Bohm and Nogales [7] for a spherical representation of particles 
still applies even to non-spherical particles providing their shape can be suitably quantified, e.g. by an 
ellipsoidal inclusion. In this particular case the Mori-Tanaka predictions were partially corroborated by 
two-dimensional numerical simulations confirming experimentally observed considerable sensitivity of 
macroscopic conductivities to the shape of particles. 

The fact that for composites with imperfect thermal contact the macroscopic predictions depend on 
particle size can be effectively handled by introducing the particle size probability density function di- 
rectly into the Mori-Tanaka estimates. Although not confirmed for material systems studied in the paper, 
this may considerably improve final predictions especially for grading curves showing significant stan- 
dard deviation of particle sizes from its mean value. This is particularly appealing, since grading curves 
are one of the few information supplied by the manufacturer. 

To conclude, it is interesting to point out that there exist many material systems that can be handled 
very effectively with simple micromechanical models with no need for laborious finite element simula- 
tions of certain representative volumes of real microstructures. 
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A. Eshelby-like tensor 

The Eshelby-like tensor for the solution of thermal conductivity problem was introduced by Hatta 
and Taya in [11]. For an ellipsoidal inclusion with semi-axes a%, a 2 , a 3 found in an isotropic matrix it 
receives the form 
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Closed form solutions of integral (39) for some special cases of ellipsoidal shapes of the inclusion can 
be found in [11]. For a general ellipsoid the solution was introduced by Chen and Yang in [25]. For 
circular and spherical shapes needed in the present study the S tensor reads 
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B. Single spherical homogeneity with imperfect interface 
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This section outlines derivation of the replacement conductivity x 1 an d the concentration factor A sp h 
introduced in Eq. (26). Since used in numerical calculations in Section 3.2, its two dimensional format 
is presented as well. It it shown that both 2D and 3D concentration factors can be recovered from the 
solution of a ID problem using a simple geometrical argument. 



Figure 11. Imperfect interface and temperature progress for: a) ID, b) 2D. 
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To that end, consider one-dimensional heat conduction problem depicted in Figure 11(a). Assum- 
ing imperfect thermal contact, the temperature drop across an infinitely thin interface layer is given by 
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Eq. (25). The local temperature gradient for perfect interface between a solitary inclusion embedded into 
an infinite matrix follows from Eq. (4) 



H 1 = A X H 



r 

x l 



H. 



(43) 



To arrive at similar format of Eq. (43) for imperfect contact, we imagine the interface temperature jump 
being smeared over the inclusion. Since the heat flux Q associated with the macroscopic temperature 
gradient H is constant throughout the composite, we obtain the total temperature change in the substitute 
inclusion in the form 
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where L stands for the inclusion length. Next, defining the local temperature gradient in the substitute 
inclusion as H 1 = A9 l / L yields the local constitutive law in terms of the replacement conductivity x 1 



Q = -x'H' 1 



-x 



L 



r 

L 



-Q 



L 2 

X 1 k 



so that 
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and in analogy with Eq. (43) we finally get 

H' 1 = A l H 
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The two-dimensional problem of a solitary circular inclusion is treated similarly. We build up on the 
assumption that the temperature gradient in the inclusion is constant and collinear with the prescribed far 
field gradient parallel to the local xi-axis, see Figure 11(b). To draw similarity with ID case we divide 
the inclusion into parallel filaments with the length L = 2a cos cp. Next, define a unit vector normal to 
the inclusion surface n = (cos cp, sin cp) T and in analogy to Eq. (44)) write the total temperature change 
in each filament for the constant heat flux q 1 = (q 1 , 0) T as 
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where d is the inclusion diameter. The equivalent conductivity has to fulfill the condition 
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Consequently, the concentration factor of the substitute inclusion attains the form 
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The analysis of a spherical inclusion follows identical steps. The replacement thermal conductivity 
for constant heat flux q 1 = (q 1 , 0, 0) T thus receives the same form as in Eq. (50) rendering the searched 
concentration factor as, recall Eqs. (27), 



A 1 = . (52) 

2x ™ + 2 
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